Inertia effects and stress accumulation in a constricted duct: 
A combined experimental and lattice Boltzmann study 
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We experimentally and numerically investigate the flow of a Newtonian fluid through a constricted 
geometry for Reynolds numbers in the range 0.1 — 100. The major aim is to study non-linear 
inertia effects at larger Reynolds numbers (>10) on the shear stress evolution in the fluid. This 
is of particular importance for blood flow as some biophysical processes in blood are sensitive to 
shear stresses, e.g., the initialization of blood clotting. We employ the lattice Boltzmann method 
for the simulations. The conclusion of the predictions is that the peak value of shear stress in the 
constriction grows disproportionally fast with the Reynolds number which leads to a non-linear shear 
stress accumulation. As a consequence, the combination of constricted blood vessel geometries and 
large Reynolds numbers may increase the risk of undesired blood clotting. 



PACS numbers: 47.ll.-j, 47.15.-x, 47.80.-v, 47.63. Cb 



I. INTRODUCTION 

There is growing evidence that blood clotting is a dy- 
namic process in which fluid shear stress plays an impor- 
tant role 0. The protein von Willebrand factor (VWF) 
shows a conformation change when the ambient shear 
rate reaches values of about 5000 s -1 . It is also known 
that arterialplaque and stents favor the emergence of 
blood clots UQ. Besides biochemical reasons, one pos- 
sible physical cause for clotting may be the detrimental 
influence of large local shear stresses in the vicinity of 
constrictions. Physically, the flow boundary conditions 
are modified by the presence of obstacles. Depending on 
the Reynolds number, obstacles can have a significant 
impact on the flow properties, even beyond the location 
of the obstacle. A prominent example is the Karman 
vortex street which is a repeating pattern of swirling vor- 
texes occurring in the laminar flow regime behind a bluff 
obstacle. 

The motivation for this article is to study inertia ef- 
fects on the shear stress in fluid flows perturbed by a 
simple obstacle and their basic implications for hemorhc- 
ology. In our present investigation, we numerically model 
and experimentally measure the flow of a Newtonian fluid 
through a duct with a simple constriction at different 
Reynolds numbers in the laminar regime (0.1 — 100). 
We emphasize that the direct comparison of experiments 
and simulations is essential. Although the accuracy and 
applicability of the employed lattice Boltzmann method 
(LBM) has been proven various times 043 ; the numer- 
ical results should be supported and verified by experi- 
ments. 
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One of the simplest symmetric, yet non-trivial, flow ge- 
ometries is a duct with a bottleneck-like constriction as 
sketched in Fig. Q] Using this obstacle, we study funda- 
mental properties of the fluid flow at different Reynolds 
numbers between 0.1 and 100. The particular design of 
the constriction has been chosen for convenience since 
this geometry is easily produced by milling techniques 
(shaping solid materials with a cutter) for the experi- 
ments. 

Since we arc interested in the basic physical effects of 
inertia in a constricted geometry, we simplify the prob- 
lem as much as possible. In particular, the fluid flow is 
assumed to be steady, and we use a Newtonian fluid (wa- 
ter in the experiments). For this reason, the Womersley 
number (a dimensionless measure for the period of pul- 
satile flows related to the viscous time scale) is zero, and 
the Reynolds number is the only relevant physical param- 
eter. The assumption of a Newtonian fluid limits the va- 
lidity of the simulations and experiments to larger blood 
vessels because at those scales the individual motion of 
the red blood cells can be neglected and the viscosity is 
virtually independent of the shear rate Q. However, this 
is no severe restriction since we are mainly interested in 
the blood flow in human coronary arteries with average 
diameters of 3-4 mm Q which is 1000 times the radius of 
a red blood cell. Typical Reynolds numbers in coronary 
arteries are of order 100 fioj ]. 

Due to the non-linear character of the Navier-Stokes 
equations, at large Reynolds numbers, abrupt changes 
in cross-section may lead to spatial variations of velocity 
and shear stress which cannot be fully understood from 
dimensional considerations or the Stokes equations. For 
this reason, we study the impact of a bottleneck-like con- 
striction on the local properties of the fluid as a function 
of the Reynolds number. We are particularly interested 
in the spatial asymmetry and the magnitude of the shear 
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Figure 1. A 2D projection of the constriction geometry used 
both in the experiments and simulations. The fluid enters 
the geometry from the left. Numerically, the inlet and outlet 
velocities u a (y,z) are taken from the analytic solution of the 
steady duct flow problem. The origin of the coordinate system 
is at the center of the constriction. The initial width (along 
y-axis) and height (along z-axis) are H, and the constricted 
width is H/2. The total length of the constriction is L C on = 
5H/4, and the initial inclination of the constriction walls is 
45° . The rounded corners with radius r = H/A are due to the 
milling technique used for fabricating the duct. The numerical 
values of the inlet and outlet duct lengths Li n and ^out are 
chosen in such a way that the flow can fully develop. 



stress. In order to emphasize the nonlinear effects arising 
at high Reynolds numbers, we vary the Reynolds number 
over three order of magnitude, (0.1 — 100), thus cover- 
ing both the fully viscous and inertial regimes. At large 
Reynolds numbers, the spatial flow velocity and shear 
stress fields are asymmetric, even in a symmetric geome- 
try. The cause is the convective term in the Navier-Stokcs 
equations. This asymmetry introduces a distinction of 
the pre- and post-constriction regions. Moreover, one 
can observe that the peak values of the shear stress close 
to the constriction increase faster than linearly with the 
Reynolds number, showing the significance of the inertia 
effects. 

The article is organized as follows. The basic hydro- 
dynamic concepts are presented in Sec. lU followed by a 
detailled description of the experimental and numerical 
setup in Sec. Mil The observations and results are pre- 
sented and discussed in Sec. HVl Finally, the conclusions 
are pointed out in Sec. [V] 



II. THEORY 

The full incompressible Navier-Stokes equations in the 
absence of a body force density read 



-Vp + 77 Am 



(1) 



where u denotes the velocity, p the density, p the pres- 
sure, and r\ the viscosity of the fluid. Introducing the 



deviatoric shear stress tensor with components 

<J a p = i] {d a Uf3 + d[)U a ) , (2) 

the viscous term in Eq. ([T]) can be written in the form 
r]Au = V ■ <t for an incompressible fluid. The devia- 
toric shear stress tensor <r (from now on only called shear 
stress) is part of the total momentum flux tensor of the 
fluid, 



M a p — p8 a p + pU a Up - (T a {}. 



(3) 



Eq. ([T|) can then be written in the compact form pdtU a = 
—df}M a fj. The first term on the right-hand-side of Eq. (J3j) 
is the isotropic pressure contribution, the second term 
denotes momentum transport due to convection (mass 
transport) which is only important at large Reynolds 
numbers. 

The shear stress a describes the momentum diffusion 
related to the viscosity of the fluid which is the central 
quantity when it comes to mechanically triggered blood 
clotting [l|. It is a second order tensor with six inde- 
pendent components (in this case only five because it is 
traceless due to the incompressibility of the fluid) . Since 
the tensor is symmetric, it has always three real eigen- 
values 03 < (T2 < <?i ■ For an incompressible fluid, the 
eigenvalues obey Tr <j = o\ + a<i + 03 = 0. In the lit- 
erature about applications of sheared fluids, usually an 
equivalent shear stress scalar a c g is provided, and the 
tensor properties are lost. Assuming that the dynamics 
of VWF is not sensitive to the tensor components cr Q( g 
but only to an effective scalar value, it arises the need for 
an appropriate definition. There arc different possibili- 
ties to construct an effective scalar measure from the full 
tensor. The von Mises stress is 



2 X! a »P a <*0 

a,/3 



(of + ff|+of) 



(4) 



<7l<7 3 . 



The last identity in Eq. ([4|) is valid if the fluid is incom- 
pressible, Tr <r = 0. The von Mises stress plays an impor- 
tant role in the analysis of materials mechanics in mate- 
rial science, but it can also be applied to fluids. Given 
the shear stress tensor at a specific point, one can also 
ask in which direction n the maximum shear component 
(T max of the tensor can be found. From a Mohr analysis 
one finds 



(5) 



The definitions in Eqs. (f4|) and §5§ and further comments 
can be found in monographs about elasticity, e.g., in 



111 ]. For an incompressible fluid, in general ovm > cr max 
holds, but one can show that if the eigenvalue 02 van- 
ishes, cr V M = fmax = o"i since then er 3 = —a\. One can 
further show that the maximum deviation between both 
shear stress scalars is about 15%. For this reason, we 
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drop a separate discussion of ct v m and <7 ma x and restrict 
ourselves to the von Mises stress. As a compact notation 
for the shear stress magnitude, we define 

a = ||er|| := a vM , (6) 

cf. Eq. ((4]). We emphasize that the concept of shear stress 
is applicable to material sciences and hydrodynamics, i.e., 
to solids and fluids. The definition for the velocity mag- 
nitude is the usual one, 



u = \\ u \\ '■= J^u a u a . (7) 

V a 

In the present article, the duct Reynolds number is 
defined as 

„ pHu , . 

Re=- . (8 

V 

H is the height and the width of the square duct, cf. Fig. 
[TJ The velocity scale u will be formally introduced in 
Eq. (|T2|) . It is the average velocity on the inlet cross- 
section with area A = H 2 and can easily be obtained in 
experiments by V = Au if the volume flux V of the fluid 
is known. 

In the limit of small Reynolds numbers and a station- 
ary situation, the left-hand-side of Eq. ([1]) is negligible, 
and one can write the Stokes equation 



Ke 


u[mms ] 


V [ml mm J 


0.1 


0.05 


0.012 


1 


0.5 


0.12 


10 


5 


1.2 


100 


50 


12 



Table I. Experimental volumetric flow rates V for achieving 
the desired Reynolds numbers and the average inlet velocities 
u. The length scale is H — 2 mm, and the fluid properties are 
p = 1000 kg m -3 and rj = 10" 3 Pas. 



Re 


r 


it 


L D /H 




0.1 


0.9 


0.000280 


0.6 


125 


1 


0.9 


0.00280 


0.6 


125 


10 


0.8 


0.0210 


0.9 


200 


20 


0.7 


0.0280 


1.4 


200 


45 


0.6 


0.0314 


2.7 


300 


100 


0.54 


0.0280 


5.8 


600 



Table II. Relevant simulation parameters: The diameter of 
the unconstricted duct in all the simulations is H — 100 lat- 
tice nodes, u is the lattice velocity at the center of the inlet 
and outlet cross-sections. It can be shown that for the cur- 
rent geometry the average velocity is u ~ 0.48m [lj. r is the 
relaxation time used for the BGK lattice Boltzmann method, 
and Ld is the development length for a given Reynolds num- 
ber according to Eq. (|16[) with D = H. Lj n and L out are the 
numerical values for the duct inlet and outlet lengths. 



= -Vp + r)Au. 



(9) 



Note that Eq. ((9|) is invariant under the transformation 
(u —> —u,\7p — > — Vp), whereas Eq. ([I]) is not. As 
a consequence, the stationary velocity field for Re = 
looks the same (up to its sign), when the flow of the 
fluid through a fixed geometry is reversed (also the signs 
of possible velocity boundary conditions have to be re- 
versed). In case of a point-symmetric geometry (invariant 
under the transformation x — > —x), as we use it here, cf. 
Fig. [IJ it is easy to see that u(x) = u(—x) must hold if 
Re = 0. The shear stress obeys <r(x) = —cr(—x) since the 
spatial derivative enters its definition, Eq. @, and leads 
to an additional minus sign. For a finite Reynolds num- 
ber, the presence of the convective term u ■ Vu breaks the 
symmetry, and the direction of the flow can be recognized 
by its inertia effects and the shape of the streamlines. 
The Karman vortex street is a demonstrative example for 
this statement: The vortexes appear only downstream. 

In this article, we study the asymmetry introduced by 
finite inertia as a function of the Reynolds number. Ad- 
ditionally, the effects of inertia on the shear stress dis- 
tribution is analyzed. For this reason, we introduce an 
index of distortion for the velocity and the shear stress. 
Knowing the analytic solutions u a (y, z ) and cr a {y 1 z) for a 
fully developed flow in the duct E2Ul4| and taking the ac- 
tual velocity Ud(y, z) and shear stress <Td(y, z) at a given 
cross-section at axial distance d from the constriction, we 



define 



y.z 



where the norm ||-|| for the shear stress and the velocity 
has been defined in Eqs. ([6]) and ([7|). The quantities 



\^2\\ u a(y,z)\\ , 



A 



y,z 



(12) 
(13) 



are the velocity and shear stress scales, defined as the av- 
erages on the inlet cross-section. They obey u oc Re and 
a oc Re and are control quantities which the numerical 
results will be related to. 

In Stokes flow, Re = 0, the fluid velocity and shear 
stress distributions look the same before and behind the 
constriction (up to signs). For this reason, Eqs. (flQ|) 
and (|lljl are invariant if the coordinate system is trans- 
formed according to x —> —x. Thus, the indexes of dis- 
tortion and I a do not change when the coordinate 
system is transformed. At finite Reynolds number, how- 
ever, the symmetry is broken, and inlet and outlet flow 
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profiles differ. Thus, we expect that I u (d) ^ I u (—d) and 
Ia(d) 7^ I a (—d) when Re ^ 0. As we will show in Section 
ITVl the slopes of I u {d) and I a {d) are well captured by a 
simple decaying exponential. For this reason, it is natu- 
ral to define a range of decay A for each exponential. It 
is given by the distance from the constriction after which 
the exponential decays to exp(— 1) of its initial value. In 
order to distinguish the ranges of decay for the velocity 
and the shear stress, we denote both quantities A„ and 
A CT , respectively. This way, the indexes of distortion can 
be approximated by 

I u {d) = I u (0)exp(-d/X u ), (14) 
I a {d) = J CT (0)exp(-d/A CT ). (15) 

Related to the definition of A u is the problem of the 
flow development length Lu in a 2D channel or a 3D 
pipe. It has been thoroughly discussed in the literature 
analytically, numerically, and experimentally due to its 
important implications in engineering [15l4l7l| . The com- 
mon approach is to impose a constant velocity profile at 
the inlet of a pipe with circular cross-section and diam- 
eter D and find the axial distance Ljj from the inlet at 
which the central velocity has reached 99% of its fully 
developed value. The parameter Ld/D is a function of 
the Reynolds number. Durst et al. [l?} have proposed 
the relation 

^ = (0.619 16 + (0.0567 Re) 16 ) 1/16 (16) 

for a Newtonian fluid in a pipe and Re = pDu/rj. This 
equation is valid for all Reynolds numbers as long as the 
flow is laminar, and the numerical error is reported to be 
< 3%. In the present simulations, there is a slightly dif- 
ferent situation: The geometry is a duct with quadratic 
cross-section, and the velocity profile at the constriction 
is not constant. However, it has turned out that Eq. (fl"6)) 
is also a good approximation for the presented problem. 
The development lengths obtained from Eq. (|16|) and us- 
ing H instead of D are shown in Tab. HU Those values 
act as a guideline for the experiments and simulations to 
assure that the duct before and behind the constriction 
is sufficiently long. 

III. EXPERIMENTAL AND NUMERICAL 
SETUP 

The employed geometry is a square duct (width = 
height = H) in the yz-plane. A 2D projection is shown in 
Fig. [I] The flow enters at the inlet in rr-direction. A con- 
striction of total length L con = 5H/A is located halfway 
between the inlet and the outlet. In the constriction, the 
width of the duct (along the y-axis) is decreased, but the 
height (along the z-axis) is not changed. The constricted 
width is H/2, leading to an average flux velocity two 
times larger than in the main duct. Due to the milling 
technique employed for the experiments, the inner edges 
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Figure 2. Locations of the velocity measurements in the ex- 
periments (black dots). All positions and distances are given 
in units of mm. The origin is located at the center of the 
constriction. The flow enters from the left. The veloci- 
ties have been measured at three positions along the :r-axis 
(x — —1.85 mm, 0.05 mm, and 1.85 mm) and five positions 
along the y-axis (y = to 0.8 mm in steps of 0.2 mm). In the 
constriction, only three data points have been taken (y — 0, 
0.2 mm, and 0.4 mm). In order to compute the velocities, the 
times of travel of the tracer particles between the dotted lines 
have been measured (Aa; = 0.8 mm and Ay = 0, respectively). 
The velocity data is shown in Tab. IHII 

of the constriction are rounded with radius r = H/4. The 
origin of the coordinate system is always located at the 
center of the constriction. 



A. Experiments 

The experimental setup consists of the duct, cf. Fig. 
[TJ a syringe pump and tubes for connecting the pump 
to the duct. The height of the duct is H = 2 mm. The 
flux is driven at a desired rate by application of a sy- 
ringe pump (NE-1000, New Era Pump Systems, Inc., NY, 
USA). Connection between the pump and the duct suc- 
ceeds over tubes connecting the syringe needle to the in- 
let. The duct itself consists of two parts. The upper part 
is made of polydimethylsiloxanc (PDMS), casted into a 
mold produced by milling. This part is then converted to 
a completely closed duct by being attached to a micro- 
scope glass slide which plays the role of the lower deck 
of the duct. Inlet and outlet are punched into the duct 
before attachment of the PDMS to the glass slide. At- 
tachment occurs by plasma oxidation of the PDMS and 
the glass slide. 

The local fluid velocity in the duct is measured by 
tracking polystyrene beads (Polysciences, Inc., War- 
rington, PA, USA) with a diameter of 10 pm which 
are suspended in the carrier fluid (water, density p = 
1000 kgnr 3 and viscosity r) = 10~ 3 Pas at 20°C). 

The experiments are conducted on a Zeiss Axiovert 
200 inverted microscope typically using a 2.5x objective. 
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The velocities of the beads are observed at a height of 
1 mm over the bottom deck of the duct (i.e., in the middle 
between bottom and top). For each Reynolds number, a 
video of the flux inside the duct is made using an ultrafast 
camera (Fastcam, Photron, CA, USA). The videos are 
analyzed by the software 'Image J' afterwards. 

Since the fluid properties and the spatial scale H are 
fixed, the Reynolds number, Eq. ([5]), can only be changed 
by choosing a mean velocity of the fluid. Four differ- 
ent Reynolds numbers have been investigated experimen- 
tally. The corresponding volume rates are shown in Tab. 

in 



B. Simulations 

The simulations were conducted with the lattice Boltz- 
mann method (LBM) using a D3Q19 BGK model [Hj]. 
There exist excellent introductory articles [l9r|2l| , mono- 
graphs 0,0, reviews 3, [22l | . and various articles about 
applications, e.g., [23l ~ [25| ] . 

In order to capture the physical boundary conditions 
of both the constriction-fluid surface (no slip) and the in- 
let and outlet cross-sections of the simulation box (fully 
developed flow), we employ the standard LBM bounce- 
back boundary condition [2l| for the former and velocity 
boundary conditions for the latter case. At the inlet and 
outlet of the computational box, a fully developed veloc- 
ity profile u a {y,z) is imposed. The velocity boundary 
condition used has been proposed by Latt et al. [l3| . We 
have chosen this approach due to its simple and straight- 
forward implementation in three-dimensional LBM sim- 
ulations. The analytic form of the stationary, fully de- 
veloped flow profile for a rectangular duct is discussed in 

ESS. 

We compute the pressure p, velocity vector u, and the 
full shear stress tensor er in the entire numerical grid. 
From this data, we can calculate the effects of inertia on 
the spatial velocity and shear stress distributions. We 
trace the maximum values of velocity and shear stress. 
Furthermore, we compute the indexes of distortion I u 
and I a , Eqs. (|10p and (fTTjl . as function of the axial dis- 
tance d from the constriction, both before and behind the 
constriction. 

In the simulations, H corresponds to 100 lattice nodes 
to ensure a sufficiently high spatial resolution. The sim- 
ulations are terminated when the relative change of ve- 
locity 




y,z) - u id(x, y,z)) (17) 

x,y,z 



becomes smaller than 10 -10 between two successive time 
steps where N is the total number of lattice nodes. This 
condition guarantees that the flow is stationary. In or- 
der to take into account the development length of a 
non-developed flow, Eq. (fT5|) . we allow the flow to relax 
towards inlet and outlet by extending the duct geometry 



correspondingly. If the inlet or outlet is too short, un- 
physical hydrodynamic interactions with the boundaries, 
such as reflections, arise. For convenience, the numerical 
inlet and outlet lengths are always identical, L m = L ont . 
The rrelevantsimulation parameters are given in Tab. [Ill 



IV. RESULTS 

a. Experimental results and comparison to simula- 
tions In Fig. [21 the locations of the velocity measure- 
ments in the experiments are presented. The velocities 
have been estimated by measuring the time of travel of 
representative tracer particles (cf. Sec. lIII Af between two 
positions along the x-axis (Aa; « 0.8 mm). The resultant 
velocity is assumed to be the velocity midway between 
the two points. The experimental data is shown together 
with the corresponding velocities from the simulations 
and the relative deviations |u eX p ~~ u s - nn \ / u s - lm in Tab. Hill 
Some streamlines found in the experiments and simula- 
tions are visualized in Figs. [4] and [5] The experimental 
and simulated results at Reynolds numbers 1 and 100 are 
shown in Fig. [4] In Fig. 03 the simulated flow fields for 
the Reynolds numbers 0.1 and 100 are directly compared. 

The experiments have been carried out with great care. 
However, as can be seen from Tab. IIII| the quantitative 
comparison of the experimental and simulation velocity 
data reveals some deviations. The major reason is that 
the velocities cannot be measured locally in the experi- 
ments. Instead, the motion of the tracer beads is followed 
over a finite distance of 0AH (0.8 mm), and the velocity 
at the middle of this line is assumed to be the average 
velocity, cf. Fig.[2[ This approach can only be accurate if 
the length over which the particles are observed is small 
compared to the typical length for the change of the ve- 
locity field. This characteristic length is of order H/4 
which is half the width of the duct inside the constric- 
tion. Consequently, there is an intrinsic uncertainty in 
the velocity measurement. The particles do not always 
move on straight lines which can be recognized from the 
shape of the streamlines, cf. Fig. [4] This makes it hard to 
achieve a good estimate for the local velocities even when 
the time resolution of the measurements is high. Espe- 
cially close to the walls (y = 0.8 mm before and behind 
and y = 0.4mm inside the constriction), the deviations 
are expected to be larger. The reason is that the veloc- 
ity gradient is maximum in the vicinity of the walls. If 
the position of the tracer beads is slightly shifted along 
the y-axis, this will lead to a large uncertainty in the 
velocity measurement. This trend can clearly be recog- 
nized in Tab. IIIII An additional, yet minor, reason for 
the deviations is that the tracer particles do not nec- 
essarily move exactly in the plane midway between the 
bottom and top walls (z = 0). However, it is encourag- 
ing to see that the qualitative shape of the experimental 
streamlines is recovered by the computer simulations. Es- 
pecially the shape of the vortexes at Re = 100, cf. Fig. 
|4(b)[ is correctly reproduced. Taking those considera- 
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position 




Re = 0.1 






Re = 1 






Re = 10 






Re = 100 




x[mm] 


\j [mm! 


exp 


sim 


dev 


exp 


sim 


dev 


exp 


sim 


dev 


exp 


sim 


dev 


— 1.85 


0.0 


0.092 


0.109 


16% 


1.2 


1.08 


11% 


9.8 


10.6 


8% 


87 


106 


18% 


-1.85 


0.2 


0.088 


0.104 


15% 


1.1 


1.03 


7% 


9.4 


10.2 


8% 


85 


101 


16% 


-1.85 


0.4 


0.076 


0.0894 


15% 


0.98 


0.891 


10% 


8.9 


8.85 


1% 


76 


88.7 


14% 


-1.85 


0.6 


0.054 


0.0669 


19% 


0.81 


0.669 


21% 


6.0 


6.69 


10% 


63 


67.1 


6% 


-1.85 


0.8 


0.038 


0.0370 


3% 


0.50 


0.371 


35% 


4.1 


3.70 


11% 


43 


35.8 


20% 


0.05 


0.0 


0.17 


0.199 


15% 


1.9 


1.99 


5% 


23 


19.4 


19% 


151 


161 


6% 


0.05 


0.2 


0.14 


0.165 


15% 


1.5 


1.65 


9% 


18 


16.1 


12% 


144 


150 


4% 


0.05 


0.4 


0.080 


0.0662 


21% 


0.93 


0.662 


40% 


6.7 


6.55 


2% 


91 


71.1 


28% 


1.85 


0.0 


0.085 


0.109 


22% 


1.0 


1.10 


9% 


14 


12.8 


9% 


144 


155 


7% 


1.85 


0.2 


0.089 


0.104 


17% 


0.92 


1.05 


12% 


10 


11.9 


16% 


140 


140 


0% 


1.85 


0.4 


0.071 


0.0894 


21% 


0.91 


0.897 


2% 


10 


9.58 


4% 


103 


95.2 


8% 


1.85 


0.6 


0.065 


0.0669 


3% 


0.59 


0.669 


12% 


5.8 


6.57 


12% 


32 


42.8 


25% 


1.85 


0.8 


0.045 


0.0370 


22% 


0.41 


0.369 


11% 


4.8 


3.35 


43% 


12 


8.46 


42% 



Table III. Measured and simulated velocities in the constriction at selected positions [x, y) midway between the bottom and 
top walls, cf. Fig. [2] All velocities are given in units of mms~ . The deviations |it e x P — Msim|/w s im are also shown. 



X u /H 



\a/H 



ric 


inlet 


outlet 


inlet 


outlet 


0.1 


0.18 


0.18 


0.20 


0.20 


1 


0.17 


0.19 


0.19 


0.21 


10 


0.16 


0.28 


0.17 


0.30 


20 


0.16 


0.44 


0.17 


0.45 


45 


0.17 


0.80 


0.17 


0.79 


100 


0.17 


1.38 


0.18 


1.37 



fies the approximations in Eqs. (Q3|) and (fT5)l and the 
introduction of the decay lengths A u and X a . 

For small Reynolds numbers (Re = 0.1 and 1), the 



Table IV. Ranges of decay X u /H and \ a jH towards the inlet 
and outlet extracted from the simulation data. The ranges 
are defined as the distances from the constriction at which 
the indexes of distortion I u and I a drop to exp( — 1) of their 
value directly at the constriction. The outlet data is also 
illustrated in Fig. [6] 



tions into account, the agreement between experiments 
and simulations is satisfactory. 

In order to test the confidence in the simulations, we 
have decreased the numerical resolution from H = 100 
to H = 80 and 40 (data not shown). We observe that 
the numerical results for H = 80 are virtually identical 
to those for H = 100 indicating that the resolution is 
sufficient to capture the correct physics. Even for H = 
40, the velocity data is accurate whereas the shear stress 
data starts to become imprecise. Due to the similarity 
of the data for H = 100 and H = 80, we believe that a 
resolution H = 100 is sufficient. In the following, we will 
only report results extracted from the simulations with 
H = 100. 

b. Index of distortion and flow relaxation In Fig. [3j 
the numerically obtained indexes of distortion for the ve- 
locity, I u , and the shear stress, I a , are presented as func- 
tion of the distance d from the constriction. The corre- 
sponding ranges of decay, X u and A CT , are shown in Tab. 
[TV] and Fig. [Sj It is obvious that the slopes of I u {d) and 
I (j (d) can be excellently described by simple exponentials 
with decay lengths X u and A CT , respectively. This justi- 



curves of I u and I a hardly depend on Re, cf. Figs. 3(a) 
and |3(b)] This is a first hint that Re = 1 still is a good 
approximation for Stokes flow. A significant change in 
the slopes is visible for larger Reynolds numbers (Re > 1) 
which can be seen from Figs. 3(c) through |3(f)| This is 
related to the influence of inertia. 

There are only small differences between the inlet and 
outlet curves of I u and I a for small Reynolds numbers, 
i.e., the flow fields are nearly symmetric with respect to 
the regions before and behind the constriction, cf. Figs. 
3(a) and |3(b)| This is another hint for the validity of the 



Stokes limit at Re < 1. For larger Re, the indexes of dis- 
tortion towards the outlet are always larger than those 
towards the inlet, indicating that the constriction mainly 
influences the flow behind itself, cf. Figs. 3(c) through 
|3(f)| Obviously, the symmetry is broken due to the pres- 
ence of inertia. This can also be seen in Fig. [7] where 
examples of the spatial shear stress evolution along the 
x-axis arc shown. For Re = 0.1, the curves are symmet- 
ric with respect to the center of the constriction, but for 
Re = 100, the asymmetry is clearly visible. 

Analyzing the data shown in Fig. [3l it is obvious that 
the decay characteristics of the distortion of the veloc- 
ity and the shear stress are similar if not identical, i.e., 
I u {d) ~ I a (d) and A u ps A^ for a given Reynolds number. 
Since the shear stress is related to the spatial derivatives 
of the velocity, this observation indicates that there is 
only one characteristic decay length both for the velocity 
and the shear stress. 

The increase of the outlet values of A u and X a with Re 
is shown in Tab. [TV] and Fig. [6] Qualitatively, the behav- 
ior of A (Re) can be understood from Eq. (fTo]) , defining 
the development length Ljj of the velocity in a pipe as 
function of the Reynolds number. Although the defini- 
tions of X u and X a on the one hand and on the other 
hand are not equivalent, both describe the same physics, 
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Figure 3. Presentation of the normalized indexes of distortion I u and I a extracted from the simulation data. The indexes of 
distortion are defined as the relative deviations of the velocity /shear stress profiles at a given cross-section with a distance d from 
the constriction compared to the reference profile at the inlet, cf. Eqs. (| 10p and |TTJ. For each Reynolds number, the normalized 
indexes of distortion are logarithmically plotted as function of the distance d/H from the beginning of the constriction to the 
inlet and from the end of the constriction to the outlet, respectively. Additionally, the exp( — l)-level is marked (horizontal, 
dotted line), defining the ranges X u /H and \ a jH . Note that the ranges on the d-axis are different for the subfigures. 



namely the relaxation behavior of the fluid as a function 
of the Reynolds number. In Stokes flow, Eq. (fTtj)) yields 
a constant development length which is also the case in 
Fig- El There is a transition region for Reynolds numbers 
in the interval [10 — 100] after which Lo(Re) becomes lin- 
ear in Rc. This linear region of \ u and X a , however, has 



not been probed in our simulations. 

From Tab. UVI we find that the ranges of decay \ u and 
\ a towards the inlet are always about 0.2H, regardless of 
the Reynolds number. The interpretation is that inertia 
affects only the fluid inside and behind the constriction. 
The fluid approaching the constriction from the inlet ex- 
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0.013 0.027 0.040 0.053 



(a)]Re = 1 and[(b) 



(b) Re = 100 

Figure 4. (Color online) The streamlines at l 
Re = 100 seen in the experiments (top) and in the simulations 
(bottom) at z = (midway between bottom and top walls), 
respectively. The fluid enters from the left. The colors in the 
simulation figures correspond to the velocity magnitudes, cf. 
Fig. \5\ 



periences the presence of the constriction only by momen- 
tum diffusion, and the Reynolds number does not play a 
significant role. This can also be seen by comparing the 
velocity and shear stress fields at different Re upstream 
of the constriction shown in Fig. [Sj For Re = 0.1 and 
100, the regions before the constriction look similar, but 
there are pronounced differences downstream. Applied 
to blood flow, this means that the constriction cannot 
cause clotting in the upstream region. 

c. Peak values of velocity and shear stress The pre- 
vious discussions clearly show that non-linear effects be- 
come important at large Reynolds numbers. In Figs. [7| 
and[8j we present additional simulation data for the ve- 
locity and the shear stress to support those observations. 




O.O^^^^T^^^g 0.00040^^053 

(a) streamlines and velocity magnitude 
0.0e+00 1.8e-05 3.6e-05 5.3e-05 7.1e-05 




0.0e^^.8e-06 3.5e-06 5.3e-06 7^-06 
(b) shear stress magnitude 

Figure 5. (Color online) Direct comparison of the simulation 
results: (a) the streamlines and velocity magnitudes and (b) 
the shear stress magnitudes at z = (midway between bottom 
and top walls) for Re = 0.1 (bottom) and Re = 100 (top), 
respectively. The fluid enters from the left. For convenience, 
the colors for the magnitudes have been chosen in such a way 
that equal values oiu/u or a /a have the same color. Numbers 
are lattice values. 
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Figure 6. The computed ranges of decay X u , A CT in the outlet 
direction from Tab. II VI are shown as function of the Reynolds 
number. For Re < 1, the ranges are constant, and the Stokes 
approximation is valid. For Re > 1, inertia effects are impor- 
tant and the relaxation of the fluid is significantly delayed. 
The results for \ u and \ a are virtually identical for a given 
Reynolds number. 
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x/H 
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(c) y/H = 0.245 

Figure 7. The shear stress a/a as & function of axial position x 
in the constriction at z = (midway between bottom and top 
walls). For the Reynolds numbers 0.1, 45, and 100, the spatial 
shear stress evolution for three different lateral positions, |(a)| 
y/H = 0.005 (close to the central axis), [(b)] y/H = 0.125 
(halfway between the central axis and the constricted wall), 
and |(c)| y/H = 0.245 (close to the constricted wall), is shown 
(y — corresponds to the central axis and y/H — 0.25 to the 
constricted wall). The center of the constriction is located at 
x = 0. The vertical dotted lines mark the beginning and end 
of the constriction. For convenience, only Re = 0.1 and 100 
are shown in |(c)| The solid vertical lines mark the averaging 
interval in order to obtain <T max /a, Eq. (|18p . For comparison, 
sec also Fig. [5(b)1 



IS 3 



a 2 







rrn r- 



- 1 1 — I — I I I I I 1 1 1 — I — I I I I 



LJJ 1_ 



0.1 



—I I I I I I I I 



Re 



10 



—i i 1111 



100 



(a) maximum velocity magnitudes 




(b) maximum shear stress magnitudes 



Figure 8. Maximum magnitudes of (a) the velocity and (b) the 
shear stress, normalized by the scales u and a defined on the 
inlet, Eqs. (|12[) and (|13|) . At small Reynolds numbers, both 
quantities do not depend on Re, and inertia is negligible. At 
larger Re, the shear stress increases while the velocity drops. 
Those are effects caused by the non-linear terms in the Navier- 
Stokes equations. 



The spatial evolution of the shear stress along the 
x-axis is shown in Fig. [7J Here, z = is fixed and 
y/H = 0.005 (close to the central axis at y = 0) in Fig. 
7(a) y/H = 0.125 (halfway between central axis and con- 

.245 (close to the 
Obviously, 



stricted walls) in Fig. 7(b) and y/H 



constricted walls at y/H = 0.25) in Fig. 7(c) 
the shear stress distribution is symmetric with respect to 
x = for Re = 0.1. At higher Reynolds numbers, a(x) 
is asymmetric. 

Another important observation is that the peak value 
of the shear stress close to the wall increases dispropor- 
tionally fast with Re, cf. Fig. 7(c) While the fluid veloc- 



ity is zero at the walls and maximum in the bulk region, 
the shear stress reaches its maximum close to or at the 
walls in a typical hydrodynamic flow situation. How- 
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ever, the lattice nature of the LBM causes inaccuracies 
in the computation of the shear stress close to inclined 
or curved obstacles. In other words: The numerical er- 
ror of the shear stress close to the wall is increased. To 
diminish this problem, we have computed the average of 
the shear stress on an interval about the maximum of the 



curve cr/cr, cf. Fig. 7(c) 



ff/8 



dx <j{x) 



(18) 



with x 2 - xi = ff/8 and (x\ + x 2 )/2 = -0.285.ff. The 
averaging process reduces possible lattice artifacts and 
is taken as a measure for the peak shear stress in the 
constriction. This procedure does not necessarily limit 
the significance of the quantity <7 max since proteins and 
cells passing the constriction do not instantaneously react 
on the local shear stress. In fact, also the time of exposure 
plays a role (which is equivalent to a finite distance along 
the path due to the advection velocity). This is well- 
known in stress induced hemolysis [26j]. The results for 
the averaged peak shear stresses tr max as a function of Re 
are presented in Fig. |8(b)| In Fig. 8(a) the maximum 
velocity u max on the centerline (y = z = 0) is shown as 
function of the Reynolds number. It is normalized by 
the characteristic velocity u to enable comparability of 
the results for different Reynolds numbers. 

In the Stokes limit, u max /u and cr max /o' do not depend 
on Re since non-linear effects are absent. The location 
where the fluid reaches its peak velocity u max is cither 
at x = (smaller Re) or behind the middle, x > 0, 
of the constriction (larger Re), cf. Fig. 5(a) At larger 
Reynolds numbers, u max /u decreases. This can be under- 
stood qualitatively by comparing the time scales for diffu- 
sion and advection. On the one hand, at small Reynolds 
numbers, a distortion in the fluid mainly propagates by 
diffusion, and advection is negligible. On the other hand, 
advection is dominant at large Reynolds numbers. Since 
the constriction is a localized perturbation at the lateral 
walls, it takes some time until it can affect the fluid in 
the vicinity of the centerline. This time can be estimated 
by the diffusion time scale 



to = 



IP 
"8V 



(19) 



where v = r\j p is the kinematic viscosity. In this time, 
however, the fluid has already propagated by a charac- 
teristic distance Ldp = uto where L^p/H oc Re. If 
Ldp is large with respect to the length of the constric- 
tion, the fluid leaves the constriction and starts to relax 
again before the fluid near the central axis is fully aware 
of the perturbation caused by the constriction. Hence, 
the centerline velocity does not as strongly increase dur- 
ing the passage through the constriction as for smaller 
Reynolds numbers. This is also the reason for the in- 
crease of <T max /<7 with Re in Fig. |8(b)| The volume flux of 
the fluid through any cross-section perpendicular to the 



x-axis has to be constant. Thus, the average fluid velocity 
must become larger inside the constriction. When the ve- 
locity near the centerline is not proportionally increased 
(which is the case at large Re) , the fluid near the walls has 
to be faster to compensate. This, on the other hand, leads 
to a disproportionate increase of the shear stress near the 
walls. In fact, for Re = 100, cr max /c7 is about 53% larger 
than in the viscous limit. This is a significant inertia ef- 
fect which will be even more severe at Re > 100. The 
implication is that unfavorable blood vessel geometries in 
combination with large Reynolds numbers can lead to a 
significant non-linear build-up of shear stress causing fur- 
ther complications during stress-induced blood clotting. 



V. CONCLUSIONS 

Large shear stresses in blood flow can lead to a confor- 
mation change of the protein von Willcbrand factor. This 
may trigger undesired blood clotting in arteries which can 
eventually lead to a coronary thrombosis. In order to es- 
timate the impact of inertia on the shear stress in coro- 
nary arteries, we have employed the lattice Boltzmann 
method to simulate the flow in a constricted geometry 
with Reynolds numbers between 0.1 and 100. We assume 
the fluid to be Newtonian since the particulate nature of 
blood and its non-Newtonian properties are only signifi- 
cant in small blood vessels like venules and arterioles. 

The major observation is that the peak value of the 
effective von Mises stress <t v m grows disproportionally 
fast with the Reynolds number in the inertial regime, 
Re > 10. At Re = 100, a common value of the Reynolds 
number in coronary arteries, the peak value of ct v m is 
more than 50% larger than expected from assuming the 
validity of Stokes flow. This observation indicates that a 
combination of pathological blood vessel geometries and 
large Reynolds numbers may increase the risk of an heart 
attack. This is a pure hydrodynamic effect. 

We further observe that the influence of the constric- 
tion is noticeable only inside and behind itself, i.e., up- 
stream of the constriction, the flow field and the shear 
stress are not significantly influenced. The downstream 
distortion decays exponentially with the distance to the 
constriction, and its range grows linear with the Reynolds 
number for large Re. In particular, the inertial effects 
break the symmetry of the flow field upstream and down- 
stream of the constriction. 

With this article, we point out that pure hydrodynamic 
effects could be the reason for an increased tendency to 
blood clotting in pathologically altered blood vessel ge- 
ometries in combination with large Reynolds numbers. 
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